%%%%%%%%%%%%%%%
% Housekeeping
%%%%%%%%%%%%%%%

clc;
clear                                   all;
close                                   all;
beep                                    off;
rng.Seed                                = 110887;
set(groot,'defaultAxesFontName', 'Garamond')
set(groot,'defaultTextFontName', 'Garamond')
set(groot,'defaultAxesTickLabelInterpreter','latex'); 
set(groot,'defaultLegendInterpreter','latex');
cd(pwd); 
addpath(strcat(pwd,'/auxiliary_functions/'));
project_folder                          = fileparts(fileparts(pwd));
root_folder                             = fileparts(project_folder );

%%%%%%%%%%%%%%%
% Model setup
%%%%%%%%%%%%%%%

% Structural parameters
pars.r                                  = -1/12*log(0.70);                                                              % Monthly model with annual discount rate of 0.70 
pars.k                                  = 0.200;                                                                        % Upper bound on adjustment intensity
pars.M_c                                = 1;                                                                            % Cash stock average
pars.M_e                                = 0.97;                                                                         % Electronic money    
pars.C                                  = 0.060;                                                                        % Strength of complementarities           
pars.mean_reversion                     = 80;                                                                           % Shock persistence, expressed as the number of days to get back to within 10% of initial value
pars.theta                              = -log(1-0.90)/( pars.mean_reversion/30 );                                      % Shock persistence, expessed as OU parameter
pars.sigma                              = 0.060;                                                                        % Standard deviation of money shock    
pars.T                                  = 12*15;                                                                        % Date after which mean-reversion goes to zero
pars.tIRF                               = 12;                                                                           % IRF horizon
S                                       = 0.20;
Alt_X0                                  = 0.10;

% Discretization parameters
Nu.NX                                   = 1e2+1;                                                                        % # of gridpoints for X
Nu.N                                    = 1e2;                                                                          % (1/2)*( # of gridpoints for M - 1 )
Nu.Nmin                                 = ceil(pars.theta/pars.sigma^2* ...
                                               ((pars.r+pars.k+pars.theta)/(pars.r+pars.k))^2* ...
                                              max((pars.M_c-pars.M_e)^2,(pars.M_e + pars.C - pars.M_c)^2)); 
Nu.N                                    = max(Nu.N,Nu.Nmin); 
Nu.Deltat                               = 1/(Nu.N*pars.theta);                                                          % Time period, expressed as fraction of a month
Nu.theta_target                         = pars.theta;
Nu.solutionmode                         = 1;

% Solve the model for different values of C
N_C                                     = 96+1;
index_baseline                          = 52;
index_noC                               = 1;
%N_C                                     = 3;
%index_baseline                          = 2;
%index_noC                               = 1;
C_vec                                   = linspace(0,0.12,N_C)';
C_vec(index_baseline)                   = 0.063;
[Compstat_a,Compstat_X]                 = deal(NaN(size(C_vec)));

for n=1:N_C

    pars.C                                 = C_vec(n);
    solution_local                         = solve(pars,Nu,1);

    if n == index_noC
        preT                                    = solution_local.preT;
        X                                       = preT.Nu.X';
        shock_index                             = find( abs(preT.Nu.M - (1-S)*preT.pars.M_c) == min(abs(preT.Nu.M - (1-S)*preT.pars.M_c)),1,'first');
        X0_index                                = find( abs(preT.Nu.X - Alt_X0) == min(abs(preT.Nu.X - Alt_X0)),1,'first');
        tvec                                    = (-preT.Nu.Deltat):preT.Nu.Deltat:preT.pars.tIRF;
        tvecIRF                                 = (-preT.Nu.Deltat):preT.Nu.Deltat:preT.pars.tIRF;
        IRF_M                                   = reshape(solution_local.IRF.IRF_M(shock_index,1,:),[solution_local.IRF.Nt+1 1]);
        IRF_X_noC                               = reshape(solution_local.IRF.IRF_X(shock_index,1,:),[solution_local.IRF.Nt+1 1]);
        IRF_a_noC                               = reshape(solution_local.IRF.IRF_a(shock_index,1,:),[solution_local.IRF.Nt+1 1]);
        Compstat_X_noC                          = reshape(solution_local.IRF.IRF_a(shock_index,:,end),[numel(X) 1]);
        Compstat_X(n)                           = IRF_X_noC(end);
        Compstat_a(n)                           = IRF_a_noC(end);
    elseif n == index_baseline
        IRF_X_C                                 = reshape(solution_local.IRF.IRF_X(shock_index,1,:),[solution_local.IRF.Nt+1 1]);
        IRF_a_C                                 = reshape(solution_local.IRF.IRF_a(shock_index,1,:),[solution_local.IRF.Nt+1 1]);
        skip                                    = 10;
        IRF_a_C_alt                             = solution_local.IRF.IRF_a(shock_index,X0_index,1:skip:end);
        Compstat_X_C                            = reshape(solution_local.IRF.IRF_a(shock_index,:,end),[numel(X) 1]);
        Compstat_X(n)                           = IRF_X_C(end);
        Compstat_a(n)                           = IRF_a_C(end);
    else
        IRF_X                                   = reshape(solution_local.IRF.IRF_X(shock_index,1,:),[solution_local.IRF.Nt+1 1]);
        IRF_a                                   = reshape(solution_local.IRF.IRF_a(shock_index,1,:),[solution_local.IRF.Nt+1 1]);
        Compstat_X(n)                           = IRF_X(end);
        Compstat_a(n)                           = IRF_a(end);
    end

end

clear                                           preT solution_local

%%%%%%%%%%
% Figure 2
%%%%%%%%%%
color1                                  = [0, 0.4470, 0.7410];
color2                                  = [0.6350, 0.0780, 0.1840];
gray                                    = [0,0,0]+0.3;
gray2                                   = [0,0,0]+0.5;
gray3                                   = [0,0,0]+0.7;
color3                                  = [0.9290, 0.6940, 0.1250];
color4                                  = [0.4660, 0.6740, 0.1880];
fignum                                  = 1;
fontsize                                = 17;

%%%%%%%%%%%%%%%
% Not exported
%%%%%%%%%%%%%%%

ff2a = figure(fignum);
fignum = fignum+1;

    % Plot
    plot(tvec,IRF_M,'-','color','black','linewidth',1); hold on;
    xlim([0 tvec(end)]);
    plot(tvec,pars.M_c*ones(size(tvec)),'--','color','black','linewidth',0.5);
    th1 = text(tvec(end)*1.05,pars.M_c,{'$M^c$'},'interpreter','latex','HorizontalAlignment','center','color','black',...
                                'VerticalAlignment','middle','fontsize',12); 
    plot(tvec(1),pars.M_c,'o','MarkerSize',5,'MarkerFaceColor','white','MarkerEdgeColor','black'); hold on;   
    plot(tvec(1),IRF_M(1),'o','MarkerSize',5,'MarkerFaceColor','black','MarkerEdgeColor','black'); hold on;
    ylim([0.8 1.2]);
    grid on;
    xlabel('$t$ (months)','interpreter','latex');
    yticks([0.8 0.9 1.0 1.1 1.2])
    title('$\mathcal{I}_{M}(t)$','interpreter','latex','fontsize',14);
    set(findall(gcf,'-property','FontSize'),'FontSize',fontsize)
 
    % Print
    set(ff2a, 'Units', 'Normalized', 'OuterPosition', [0, 0, 0.3, 0.3725]);    
    set(ff2a,'Units','inches');
    screenposition                              = get(ff2a,'Position');
    set(ff2a,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2a.PaperPositionMode                       ='auto';
    %exportgraphics(ff2a,[root_folder '/output/figures/' 'IRF_M.pdf'],'BackgroundColor','none');  

%%%%%%%%%%%%
% Figure 4-2
%%%%%%%%%%%%

ff2b = figure(fignum);
fignum = fignum+1;

    % Plot
    plot(C_vec,Compstat_X,'-x','Color','black','linewidth',1); hold on;
    xlim([0 max(C_vec)]);
    ylim([0 1]);
    grid on;
    hold on;
    xlabel('$C$','interpreter','latex');
    title('$\mathcal{I}_{X}(t;0,C)$','interpreter','latex');
    set(findall(gcf,'-property','FontSize'),'FontSize',fontsize)
    
    % Print
    set(ff2b, 'Units', 'Normalized', 'OuterPosition', [0, 0, 0.3, 0.3725]);    
    set(ff2b,'Units','inches');
    screenposition                              = get(ff2b,'Position');
    set(ff2b,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2b.PaperPositionMode                       ='auto';
    exportgraphics(ff2b,[root_folder '/output/figures/' 'Figure_4_2.pdf'],'BackgroundColor','none');  

%%%%%%%%%%%%
% Figure 4-1
%%%%%%%%%%%%

ff2c = figure(fignum);
fignum = fignum+1;

    % Plot
    p1=plot(tvec,IRF_X_C,'-','color',color2,'linewidth',1); hold on;
    p2=plot(tvec,IRF_X_noC,'-','color',color1,'linewidth',1); 
    xlim([0 tvec(end)]);
    ylim([0 0.4]);
    grid on;
    plot(tvec(1),IRF_X_C(1),'o','MarkerSize',5,'MarkerFaceColor',color1,'MarkerEdgeColor',color1); hold on;
    xlabel('$t$ (months)','interpreter','latex');
    legend([p1,p2],{'$C > 0$','$C = 0$'},'interpreter','latex','location','Northwest');
    legend boxoff;
    title('$\mathcal{I}_{X}(t;0;C)$','interpreter','latex','fontsize',14);
    set(findall(gcf,'-property','FontSize'),'FontSize',fontsize)

    % Print
    set(ff2c, 'Units', 'Normalized', 'OuterPosition', [0, 0, 0.3, 0.3725]);    
    set(ff2c,'Units','inches');
    screenposition                              = get(ff2c,'Position');
    set(ff2c,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2c.PaperPositionMode                       ='auto';
    exportgraphics(ff2c,[root_folder '/output/figures/' 'Figure_4_1.pdf'],'BackgroundColor','none');

%%%%%%%%%%%%
% Figure 4-4
%%%%%%%%%%%%

ff2d = figure(fignum);
fignum = fignum+1;

    % Plot
    plot(C_vec,Compstat_a,'-x','Color','black','linewidth',1); hold on;
    xlim([0 max(C_vec)]);
    ylim([0 1]);
    grid on;
    hold on;
    xlabel('$C$','interpreter','latex');
    title('$\mathcal{I}_{a}(t;0,C)$','interpreter','latex');
    set(findall(gcf,'-property','FontSize'),'FontSize',fontsize)
    
    % Print
    set(ff2d, 'Units', 'Normalized', 'OuterPosition', [0, 0, 0.3, 0.3725]);    
    set(ff2d,'Units','inches');
    screenposition                              = get(ff2d,'Position');
    set(ff2d,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2d.PaperPositionMode                       ='auto';
    exportgraphics(ff2d,[root_folder '/output/figures/' 'Figure_4_4.pdf'],'BackgroundColor','none');  
    
%%%%%%%%%%%%
% Figure 4-3
%%%%%%%%%%%%

ff2e = figure(fignum);
fignum = fignum+1;

    % Plot
    p1=plot(tvec,IRF_a_C,'-','color',color2,'linewidth',1); hold on;
    p2=plot(tvec,IRF_a_noC,'-','color',color1,'linewidth',1); 
    xlim([0 tvec(end)]);
    ylim([0 1]);
    grid on;
    plot(tvec(1),IRF_a_C(1),'o','MarkerSize',5,'MarkerFaceColor',color1,'MarkerEdgeColor',color1); hold on;
    xlabel('$t$ (months)','interpreter','latex');
    legend([p1,p2],{'$C > 0$','$C = 0$'},'interpreter','latex','location','Northeast');
    legend boxoff;
    title('$\mathcal{I}_{a}(t;0,C)$','interpreter','latex','fontsize',14);
    set(findall(gcf,'-property','FontSize'),'FontSize',fontsize)

    % Print
    set(ff2e, 'Units', 'Normalized', 'OuterPosition', [0, 0, 0.3, 0.3725]);    
    set(ff2e,'Units','inches');
    screenposition                              = get(ff2e,'Position');
    set(ff2e,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2e.PaperPositionMode                       ='auto';
    exportgraphics(ff2e,[root_folder '/output/figures/' 'Figure_4_3.pdf'],'BackgroundColor','none'); 

%%%%%%%%%%%%
% Figure 4-6
%%%%%%%%%%%%

ff2f = figure(fignum);
fignum = fignum+1;

    p1=plot(X,Compstat_X_C,'-','color',color2,'linewidth',1); hold on;
    p2=plot(X,Compstat_X_noC,'-','color',color1,'linewidth',1); hold on;
    xlim([0 max(X)]);
    ylim([0 1]);
    grid on;
    hold on;
    xlabel('$X_0$','interpreter','latex');
    title({'$\mathcal{I}_{a}(t;X_0,C)$'},'interpreter','latex');
    legend([p1,p2],{'$C > 0$','$C = 0$'},'interpreter','latex','location','East');
    legend boxoff;
    set(findall(gcf,'-property','FontSize'),'FontSize',fontsize)

    % Print
    set(ff2f, 'Units', 'Normalized', 'OuterPosition', [0, 0, 0.3, 0.3725]);    
    set(ff2f,'Units','inches');
    screenposition                              = get(ff2f,'Position');
    set(ff2f,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2f.PaperPositionMode                       ='auto';
    exportgraphics(ff2f,[root_folder '/output/figures/' 'Figure_4_6.pdf'],'BackgroundColor','none'); 

%%%%%%%%%%%%
% Figure 4-5
%%%%%%%%%%%%

ff2g = figure(fignum);
fignum = fignum+1;

    % Plot
    p1=plot(tvec,IRF_a_C,'-','color',color2,'linewidth',1); hold on;
    p2=plot(tvec(1:skip:end),reshape(IRF_a_C_alt,[numel(IRF_a_C_alt) 1]),'-x','color',color2,'linewidth',1); hold on; 
    xlim([0 tvec(end)]);
    ylim([0 1]);
    grid on;
    plot(tvec(1),IRF_a_C(1),'o','MarkerSize',5,'MarkerFaceColor',color1,'MarkerEdgeColor',color1); hold on;
    xlabel('$t$ (months)','interpreter','latex');
    legend([p1,p2],{'$C > 0$, $X_0=0$','$C > 0$, $X_0=0.1$'},'interpreter','latex','location','Northeast');
    legend boxoff;
    title('$\mathcal{I}_{a}(t;X_0,C)$','interpreter','latex','fontsize',14);
    set(findall(gcf,'-property','FontSize'),'FontSize',fontsize)

    % Print
    set(ff2g, 'Units', 'Normalized', 'OuterPosition', [0, 0, 0.3, 0.3725]);    
    set(ff2g,'Units','inches');
    screenposition                              = get(ff2g,'Position');
    set(ff2g,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2g.PaperPositionMode                       ='auto';
    exportgraphics(ff2g,[root_folder '/output/figures/' 'Figure_4_5.pdf'],'BackgroundColor','none');     